---
title: "Step4: regression analysis"
author: "Endre Borbáth and Swen Hutter"
output: html_document
---

# Loading packages

```{r, echo=T, message = T, include = T}

library(readxl)
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(grid)
library(haven)
library(forcats)
library(lmtest)
library(plm)
library(texreg)
library(lme4)
library(broom.mixed)
library(performance)
library(kableExtra)
library(thematic)
library(effects)
library(fixest)
library(parameters)
library(margins)
library(here)

rm(list = ls()) 


load("dataset_long.Rdata")

range01 <- function(x){(x-min(x, na.rm = TRUE))/(max(x, na.rm = TRUE)-min(x, na.rm = TRUE))}

colors <- thematic::okabe_ito(8)[-6]

# to calculate margins with 90 and 95 percent CIs
margins_2ci <- function(mod) { 

mar <- margins(mod, type="response")

plot_dat <- summary(mar, level = 0.95) %>% 
  rename(`2.5 %` = lower, 
         `97.5 %` = upper) 

plot_dat2 <- summary(mar, level = 0.90) %>% 
  select(lower, upper) %>% 
  rename(`5 %` = lower,
         `95 %` = upper)

plot_dat <- bind_cols(plot_dat, plot_dat2)
  
plot_dat$model <- as.character(deparse(substitute(mod)))

return(plot_dat)

}


```

# Appendix A

## Figure 7: Coverage of the VParty variables

```{r, echo=T, message = T, include = T}

# Step 1: Identify which (country, year) pairs exist in the original data
existing_cells <- d_lng %>%
  distinct(country, year)

# Step 2: Determine if each existing cell is fully NA or not
cell_status <- d_lng %>%
  group_by(country, year) %>%
  summarize(all_na = all(is.na(VDEM_v2panom_lag1)), .groups = "drop") %>%
  mutate(status = ifelse(all_na, "expected_but_missing", "observed"))

# Step 3: Combine with all possible (country, year) combos to identify missing-by-design
full_grid <- expand_grid(
  country = sort(unique(d_lng$country)),
  year = sort(unique(d_lng$year))
) %>%
  inner_join(existing_cells, by = c("country", "year")) %>%  # only keep expected cells
  left_join(cell_status, by = c("country", "year"))

# Fill in NA (i.e., expected but no data) → should be "expected_but_missing"
full_grid <- full_grid %>%
  mutate(status = ifelse(is.na(status), "expected_but_missing", status))

# Step 4: Plot with fill + outline logic
ggplot(full_grid, aes(x = year, y = country)) +
  geom_tile(data = full_grid %>% filter(status == "observed"),
            fill = "gray50", color = NA) +
  geom_tile(data = full_grid %>% filter(status == "expected_but_missing"),
            fill = NA, color = "black") +
  scale_x_continuous(
    breaks = sort(unique(d_lng$year)),
    expand = c(0, 0)
  ) +
  scale_y_discrete(limits = rev) +
  labs(x = NULL, y = NULL) +
  theme_linedraw(base_size = 11) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    panel.grid = element_blank(),
    legend.position = "none"
  )

ggsave(plot=last_plot(),
       filename="appendix/figures/missing_map.png",
       width=14.5,
       height=14.5,
       scale=1.4,
       dpi=400,
       units="cm")

```

# Full sample

- In this part I don't yet add the party organizational variables, since those have more missings

```{r, echo=T, message = T, include = T}

party_lev_dummy <- c("past_pm_lag1", "opp_party_lag1", "coalition")
party_lev_cont <- c("vote_share_lag1", "party_age") # "state_market", "liberty_authority"
party_lev_with_missings <- c("VDEM_v2panom_lag1", "VDEM_v2padisa_lag1")
context_lev <- c("cmp_econ_polar_lag1",
                 "cmp_cult_polar_lag1", "tv_it_lag1", "Voter_Turnout_lag1",
                 "v2x_polyarchy_lag1")
context_lev_dummy <- "gtm_parl_lag1"
dep_vars <- c("dep_var1", "dep_var2", "dep_var3")
id_vars <- c("party_id", "year", "country_year", "country")

# party and party year variables

reg_dat <- d_lng %>% 
  # filter(democratization=="Ist or IInd wave") %>%
  mutate(dep_var1 = case_when(typology == "Ideological party" ~ 1,
                              typology != "Ideological party" ~ 0),
         dep_var2 = case_when(typology == "Nonideological nonparty" ~ 1,
                                  typology != "Nonideological nonparty" ~ 0),
         dep_var3 = case_when(typology == "Nonideological movement" ~ 1,
                                  typology != "Nonideological movement" ~ 0),
         party_fam = case_when(family_name %in% c("to be coded", 
                                                  "Special issue", 
                                                  "no family") ~ "Other",
                               TRUE ~ family_name),
         opp_party_lag1 = case_when(cabinet_party_lag1 == 0 ~ 1,
                                    cabinet_party_lag1 == 1 ~ 0),
         party_age = log(party_age), # take the log before the centering
         gtm_parl_lag1 = factor(gtm_parl_lag1, levels=c(0, 1)),
         country_year=paste(country, year, sep="_")) %>% 
  mutate(decades=case_when(year<1970 ~ "1945-1969",
                           year>=1970 & year<1990 ~ "1970-1989",
                           year>=1990 & year<2009 ~ "1990-2008",
                           year>=2009 & year<=2023 ~ "2009-2023")) %>% 
  mutate(decades=factor(decades, levels = c("1945-1969", "1970-1989",
                                            "1990-2008", "2009-2023"))) %>% 
  mutate(party_fam=factor(party_fam, levels=c("Christian democracy",  
                                              "Liberal", 
                                              "Social democracy",
                                              "Communist/Socialist", 
                                              "Green/Ecologist", 
                                              "Conservative", 
                                              "Right-wing",
                                              "Agrarian",
                                              "Other"))) %>% 
  mutate(country_decades=paste0(country, "_", as.character(decades))) %>% 
  # mutate_at(vars(all_of(party_lev_dummy)), ~ factor(., levels = c("0", "1"))) %>% 
  mutate_at(vars(coalition, opp_party_lag1), ~ factor(., levels = c("0", "1"))) %>% 
  mutate_at(vars(all_of(c(party_lev_cont, party_lev_with_missings))), ~ range01(as.numeric(.)))

# country, country year and year level variables

reg_dat <- reg_dat %>%
  mutate_at(vars(all_of(context_lev), year), ~range01(as.numeric(.))) 

# centering, and selecting what is included
reg_dat <- reg_dat %>%  
  select(all_of(party_lev_dummy), all_of(party_lev_cont), all_of(party_lev_with_missings),
         all_of(context_lev), all_of(dep_vars), all_of(id_vars), 
         decades, year, party_fam, coalition, country_decades,
         all_of(context_lev_dummy), typology) %>% 
  group_by(country_year) %>%
  mutate_at(vars(all_of(c(party_lev_cont, party_lev_with_missings))), ~.-mean(., na.rm=TRUE)) %>%
  ungroup(.) %>%
  group_by(country) %>% 
  mutate_at(vars(all_of(c(context_lev))), ~.-mean(unique(., na.rm=TRUE), na.rm=TRUE)) %>% 
  ungroup(.) %>% 
  # mutate_at(vars(year), ~.-mean(unique(., na.rm=TRUE), na.rm=TRUE)) %>% 
  drop_na(all_of(party_lev_dummy), all_of(party_lev_cont), all_of(context_lev),
          all_of(context_lev_dummy), all_of(dep_vars), all_of(id_vars))

#write_dta(reg_dat, "reg_dat.dta")

```


```{r, echo=T, message = T, include = T}


m.full.1 <- glmer(dep_var1 ~
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 + 
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

m.full.2 <- glmer(dep_var2 ~ 
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 + 
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

m.full.3 <- glmer(dep_var3 ~ 
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 +
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

```

# Appendix A

## Table 2: Hierarchical logistic regression models of party brands (1945-2023)

```{r, echo=T, message = T, include = T}

latex_table <- texreg(list(m.full.1, m.full.2, m.full.3),
       custom.model.names = c("party name", 
                              "nonparty name", 
                              "movement name"),
       custom.header = list("Classical" = 1, "Nonideological" = 2, "Nonideological"=3),
       doctype = FALSE, center = TRUE, dcolumn=TRUE, digits=2,
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!", return.string = TRUE,
       single.row=TRUE, no.margin = TRUE,
       caption = "Hierarchical logistic regression models of party brands (1945-2023)",
       caption.above = TRUE,
       label = "table:mlm_reg_table_all",
       include.loglik = FALSE, include.variance = FALSE,
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "party_famCommunist/Socialist" = "Communist/Socialist",
                              "party_famSocial democracy" = "Social democracy",
                              "party_famGreen/Ecologist" = "Green/Ecologist",
                              "party_famLiberal" = "Liberal",
                              "party_famConservative" = "Conservative",
                              "party_famAgrarian" = "Agrarian",
                              "party_famRight-wing" = "Right-wing",
                              "party_famOther" = "Other",
                              "past_pm_lag1has been/is pm party" = "Past PM party",
                              "opp_party_lag11" = "Opp. party",
                              "party_age" = "Log(party age)",
                              "vote_share_lag1" = "Vote share",
                              "coalition1" = "Electoral coalition",
                              "cmp_econ_polar_lag1" = "Econ. polariz.",
                              "cmp_cult_polar_lag1" = "Cult. polariz.",
                              "tv_it_lag1" = "Electoral volatil.",
                              "Voter_Turnout_lag1" = "Turnout",
                              "v2x_polyarchy_lag1" = "Electoral dem.",
                              "gtm_parl_lag11" = "Presidentialism",
                              "decades1970-1989" = "1970-1989",
                              "decades1990-2008" = "1990-2008",
                              "decades2009-2023" = "2009-2023"),
       groups = list("Party family (ref: Christ. dem.)"=2:9,
                     "Gov. status" = 10:11,
                     "Party features"=12:14,
                     "System features"=15:20,
                     "Decades (ref: 1945-1969)"=21:23),
       custom.gof.names	= c(NA, NA, NA, "Country*Year", "Country"))
  
# Remove the offending \midrule.
# (Adjust the pattern if needed to precisely target the unwanted line.)
latex_table <- gsub("(?m)^\\\\cmidrule.*\n", "", latex_table, perl = TRUE)

# Write the modified LaTeX code to file.
writeLines(latex_table, con = "appendix/tables/mlm_logits.tex")


```


# Main text

## Figure 4: Logistic regression model of party naming

```{r, echo=T, message = T, include = T}

plot_dat <- rbind(margins_2ci(m.full.1), 
                  margins_2ci(m.full.2), 
                  margins_2ci(m.full.3))

plot_dat <- plot_dat %>% 
  rename(estimate=AME,
         term=factor) %>% 
  filter(!(grepl("(Intercept)", term))) %>% 
  mutate(DV=case_when(model %in% c("m.full.1") ~ "Classical\nparty name",
                      model %in% c("m.full.2") ~ "Nonideological\nnonparty name",
                      model %in% c("m.full.3") ~ "Nonideological\nmovement name")) %>% 
  mutate(DV=factor(DV, levels = c("Classical\nparty name", 
                              "Nonideological\nnonparty name", 
                              "Nonideological\nmovement name"))) %>% 
  mutate(DV=fct_rev(DV)) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(term=case_when(
    term=="party_famAgrarian" ~ "Agrarian",
    term=="party_famCommunist/Socialist" ~ "Communist/Socialist",
    term=="party_famConservative" ~ "Conservative",
    term=="party_famGreen/Ecologist" ~ "Green/Ecologist",
    term=="party_famLiberal" ~ "Liberal",
    term=="party_famOther" ~ "Other",
    term=="party_famRight-wing" ~ "Right-wing",
    term=="party_famSocial democracy" ~ "Social democracy", 
    term=="past_pm_lag1has been/is pm party" ~ "Past PM party",
    term=="opp_party_lag11" ~ "Opp. party",
    term=="party_age" ~ "Log(party age)",
    term=="vote_share_lag1" ~ "Vote share",
    term=="coalition1" ~ "Electoral coalition",
    term=="Voter_Turnout_lag1" ~ "Turnout",
    term=="cmp_econ_polar_lag1" ~ "Econ. polariz.",
    term=="cmp_cult_polar_lag1" ~ "Cult. polariz.",
    term=="tv_it_lag1" ~ "Electoral volatility",
    term=="v2x_polyarchy_lag1" ~ "Electoral dem.",
    term== "gtm_parl_lag11" ~ "Presidentialism",
    term=="decades1970-1989" ~ "1970-1989",
    term=="decades1990-2008" ~ "1990-2008",
    term=="decades2009-2023" ~ "2009-2023"
  )) %>% 
  mutate(term=factor(term, levels = c("Communist/Socialist", "Social democracy",
                                      "Green/Ecologist", "Liberal",
                                      "Conservative", "Agrarian", 
                                      "Right-wing", "Other", 
                                      "Past PM party", "Log(party age)", "Opp. party",
                                      "Vote share", "Electoral coalition", 
                                      "Econ. polariz.", "Cult. polariz.",
                                      "Electoral volatility", "Turnout", 
                                      "Electoral dem.", "Presidentialism",
                                      "1970-1989", "1990-2008", "2009-2023"))) %>% 
  mutate(term=fct_rev(term)) %>%
  # mutate_at(vars(estimate, `2.5 %`, `97.5 %`, `5 %`, `95 %`), ~exp(.)) %>%
  mutate(facets_var=case_when(term %in% c("Liberal", 
                                      "Social democracy", "Communist/Socialist", 
                                      "Green/Ecologist", "Conservative", 
                                      "Right-wing", "Agrarian", 
                                      "Other") ~ "Party family (ref: Christian democracy)", 
                          term %in% c("Past PM party", "Opp. party", "Log(party age)", 
                                      "Vote share", "Electoral coalition") ~ "Party features",
                          term %in% c("Turnout",  "Econ. polariz.", 
                                      "Cult. polariz.", "Electoral volatility",
                                      "Electoral dem.", "Presidentialism") ~ "System features",
                          term %in% c("1970-1989", "1990-2008", 
                                      "2009-2023") ~ "Decades (ref: 1945-1969)")) %>% 
  mutate(facets_var=factor(facets_var, levels = c("Party family (ref: Christian democracy)",
                                          "Party features", "System features",
                                          "Decades (ref: 1945-1969)"))) %>% 
  select(DV, term, estimate, facets_var, `2.5 %`, `97.5 %`, `5 %`, `95 %`) 

GeomPointrange$draw_key <-  function (data, params, size)     {
  
  draw_key_vpath <- function (data, params, size) {
    segmentsGrob(0.1, 0.5, 0.9, 0.5, 
                 gp = gpar(col = alpha(data$colour, data$alpha), 
                           lwd = data$size * .pt, lty = data$linetype, 
                           lineend = "butt"), arrow = params$arrow)
  }
  
  grid::grobTree(draw_key_vpath(data, params, size), 
                 draw_key_point(transform(data, size = data$size * 4), params))
}


space_between_bars <- 0.7

ggplot(plot_dat, aes(y=term, color=DV)) +
  geom_vline(xintercept = 0, linetype="dashed", color="gray60") +
  geom_point(aes(x=estimate, shape=DV), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes(xmin=`2.5 %`, xmax=`97.5 %`,
  group=DV), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=`5 %`, xmax=`95 %`,
  group=DV), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("Average Marginal Effect (Response Scale)") +
  facet_wrap(~facets_var, ncol=2, scales = "free_y") + #
  scale_color_manual(name="", values=c("#0072B2", "#56B4E9", "#D55E00"), 
                     guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 15)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(2,"cm"),
        axis.title.y = element_blank(),
        axis.text.y = element_text(hjust=0.5))

ggsave(plot=last_plot(),
       filename="main_text/figures/Figure4.png",
       width=14.5,
       height=14.5,
       scale=1.2,
       dpi=400,
       units="cm")


```


# Partial sample

```{r, echo=T, message = T, include = T}

# centering, and selecting what is included
reg_dat <- reg_dat %>%  
  drop_na(all_of(c(party_lev_with_missings)))

```


```{r, echo=T, message = T, include = T}

m.partial.1 <- glmer(dep_var1 ~
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 +
                    VDEM_v2panom_lag1 + VDEM_v2padisa_lag1 +
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

m.partial.2 <- glmer(dep_var2 ~ 
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 + 
                    VDEM_v2panom_lag1 + VDEM_v2padisa_lag1 +
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

m.partial.3 <- glmer(dep_var3 ~ 
                    party_fam + past_pm_lag1 + opp_party_lag1 + 
                    vote_share_lag1 + party_age + 1 + 
                    cmp_econ_polar_lag1 + cmp_cult_polar_lag1 + 
                    VDEM_v2panom_lag1 + VDEM_v2padisa_lag1 +
                    Voter_Turnout_lag1 + tv_it_lag1 + 
                    v2x_polyarchy_lag1 + gtm_parl_lag1 +
                    coalition + decades +
                    (1 | country_year) + (1 | country),
                  family = binomial(link='logit'),
                  data=reg_dat)

```


# Appendix A

## Table 3: Partial sample (1970-2023) - hierarchical logistic regression models of party brands with additional independent variables

```{r, echo=T, message = T, include = T}

latex_table <- texreg(list(m.partial.1, m.partial.2, m.partial.3),
       custom.model.names = c("party name", 
                              "nonparty name", 
                              "movement name"),
       custom.header = list("Classical" = 1, "Nonideological" = 2, "Nonideological"=3),
       doctype = FALSE, center = TRUE, dcolumn=TRUE, digits=2,
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!", return.string = TRUE,
       single.row=TRUE, no.margin = TRUE,
       caption.above = TRUE,
       caption = "Partial sample (1970-2023) - hierarchical logistic regression models of party brands with additional independent variables",
       label = "table:mlm_reg_table_partial",
       include.loglik = FALSE, include.variance = FALSE,
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "party_famRight-wing" = "Right-wing",
                              "party_famGreen/Ecologist" = "Green/Ecologist",
                              "party_famCommunist/Socialist" = "Communist/Socialist",
                              "party_famConservative" = "Conservative",
                              "party_famSocial democracy" = "Social democracy",
                              "party_famLiberal" = "Liberal",
                              "party_famAgrarian" = "Agrarian",
                              "party_famOther" = "Other",
                              "VDEM_v2panom_lag1" = "Candidate nomination",
                              "VDEM_v2padisa_lag1" = "Internal cohesion",
                              "past_pm_lag1has been/is pm party"= "Past PM party",
                              "opp_party_lag11" = "Opp. party",
                              "party_age" = "Log(party age)",
                              "vote_share_lag1" = "Vote share",
                              "coalition1" = "Electoral coalition",
                              "cmp_econ_polar_lag1" = "Econ. polariz.",
                              "cmp_cult_polar_lag1" = "Cult. polariz.",
                              "tv_it_lag1" = "Electoral volatil.",
                              "Voter_Turnout_lag1" = "Turnout",
                              "v2x_polyarchy_lag1" = "Electoral dem.",
                              "gtm_parl_lag11" = "Presidentialism",
                              "decades1990-2008" = "1990-2008",
                              "decades2009-2023" = "2009-2023"),
       groups = list("Party family (ref: Christ. dem.)"=2:9,
                     "Party features"=10:16,
                     "System features"=17:22,
                     "Decades (ref: 1970-1989)"=23:24),
       custom.gof.names	= c(NA, NA, NA, "Country*Year", "Country"))

# Remove the offending \midrule.

# (Adjust the pattern if needed to precisely target the unwanted line.)
latex_table <- gsub("(?m)^\\\\cmidrule.*\n", "", latex_table, perl = TRUE)

# Write the modified LaTeX code to file.
writeLines(latex_table, con = "appendix/tables/mlm_logits2.tex")

```

## Figure 8: Logistic regression model of party naming

```{r, echo=T, message = T, include = T}

plot_dat <- rbind(margins_2ci(m.partial.1), 
                  margins_2ci(m.partial.2), 
                  margins_2ci(m.partial.3))


plot_dat <- plot_dat %>% 
  rename(estimate=AME,
         term=factor) %>% 
  filter(!(grepl("(Intercept)", term))) %>% 
  mutate(DV=case_when(model %in% c("m.partial.1") ~ "Classical\nparty name",
                      model %in% c("m.partial.2") ~ "Nonideological\nnonparty name",
                      model %in% c("m.partial.3") ~ "Nonideological\nmovement name")) %>% 
  mutate(DV=factor(DV, levels = c("Classical\nparty name", 
                              "Nonideological\nnonparty name", 
                              "Nonideological\nmovement name"))) %>% 
  mutate(DV=fct_rev(DV)) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(term=case_when(
    term=="party_famAgrarian" ~ "Agrarian",
    term=="party_famCommunist/Socialist" ~ "Communist/Socialist",
    term=="party_famConservative" ~ "Conservative",
    term=="party_famGreen/Ecologist" ~ "Green/Ecologist",
    term=="party_famLiberal" ~ "Liberal",
    term=="party_famOther" ~ "Other",
    term=="party_famRight-wing" ~ "Right-wing",
    term=="party_famSocial democracy" ~ "Social democracy", 
    term=="VDEM_v2panom_lag1" ~ "Candidate nomination",
    term=="VDEM_v2padisa_lag1" ~ "Internal cohesion",
    term=="past_pm_lag1has been/is pm party" ~ "Past PM party",
    term=="opp_party_lag11" ~ "Opp. party",
    term=="party_age" ~ "Log(party age)",
    term=="vote_share_lag1" ~ "Vote share",
    term=="coalition1" ~ "Electoral coalition",
    term=="Voter_Turnout_lag1" ~ "Turnout",
    term=="cmp_econ_polar_lag1" ~ "Econ. polariz.",
    term=="cmp_cult_polar_lag1" ~ "Cult. polariz.",
    term=="tv_it_lag1" ~ "Electoral volatility",
    term=="v2x_polyarchy_lag1" ~ "Electoral dem.",
    term== "gtm_parl_lag11" ~ "Presidentialism",
    term=="decades1970-1989" ~ "1970-1989",
    term=="decades1990-2008" ~ "1990-2008",
    term=="decades2009-2023" ~ "2009-2023"
  )) %>% 
  mutate(term=factor(term, levels = c("Communist/Socialist", "Social democracy",
                                      "Green/Ecologist", "Liberal",
                                      "Conservative", "Agrarian", 
                                      "Right-wing", "Other", 
                                      "Candidate nomination", "Internal cohesion",
                                      "Past PM party", "Log(party age)", "Opp. party",
                                      "Vote share", "Electoral coalition", 
                                      "Econ. polariz.", "Cult. polariz.",
                                      "Electoral volatility", "Turnout", 
                                      "Electoral dem.", "Presidentialism",
                                      "1970-1989", "1990-2008", "2009-2023"))) %>% 
  mutate(term=fct_rev(term)) %>%
  # mutate_at(vars(estimate, `2.5 %`, `97.5 %`, `5 %`, `95 %`), ~exp(.)) %>%
  mutate(facets_var=case_when(term %in% c("Liberal", 
                                      "Social democracy", "Communist/Socialist", 
                                      "Green/Ecologist", "Conservative", 
                                      "Right-wing", "Agrarian", 
                                      "Other") ~ "Party family (ref: Christian democracy)", 
                          term %in% c("Candidate nomination", "Internal cohesion",
                                      "Past PM party", "Opp. party", "Log(party age)", 
                                      "Vote share", "Electoral coalition") ~ "Party features",
                          term %in% c("Turnout",  "Econ. polariz.", 
                                      "Cult. polariz.", "Electoral volatility",
                                      "Electoral dem.", "Presidentialism") ~ "System features",
                          term %in% c("1970-1989", "1990-2008", 
                                      "2009-2023") ~ "Decades (ref: 1945-1969)")) %>% 
  mutate(facets_var=factor(facets_var, levels = c("Party family (ref: Christian democracy)",
                                          "Party features", "System features",
                                          "Decades (ref: 1945-1969)"))) %>% 
  select(DV, term, estimate, facets_var, `2.5 %`, `97.5 %`, `5 %`, `95 %`) 

GeomPointrange$draw_key <-  function (data, params, size)     {
  
  draw_key_vpath <- function (data, params, size) {
    segmentsGrob(0.1, 0.5, 0.9, 0.5, 
                 gp = gpar(col = alpha(data$colour, data$alpha), 
                           lwd = data$size * .pt, lty = data$linetype, 
                           lineend = "butt"), arrow = params$arrow)
  }
  
  grid::grobTree(draw_key_vpath(data, params, size), 
                 draw_key_point(transform(data, size = data$size * 4), params))
}


space_between_bars <- 0.7

ggplot(plot_dat, aes(y=term, color=DV)) +
  geom_vline(xintercept = 0, linetype="dashed", color="gray60") +
  geom_point(aes(x=estimate, shape=DV), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes(xmin=`2.5 %`, xmax=`97.5 %`,
  group=DV), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=`5 %`, xmax=`95 %`,
  group=DV), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("Average Marginal Effect (Response Scale)") +
  facet_wrap(~facets_var, ncol=2, scales = "free_y") + #
  scale_color_manual(name="", values=c("#0072B2", "#56B4E9", "#D55E00"), 
                     guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 15)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(2,"cm"),
        axis.title.y = element_blank(),
        axis.text.y = element_text(hjust=0.5))

ggsave(plot=last_plot(),
       filename="appendix/figures/logit_partial_sample.png",
       width=14.5,
       height=14.5,
       scale=1.4,
       dpi=400,
       units="cm")

```


```{r, echo=T, message = T, include = T}

```